This R markdown file contains the code to reproduce the figures in “Hydraulic prediction of drought-induced plant dieback and top-kill depends on leaf habit and growth form” by Chen Ya-Jun et al. in Ecology Letters.

1 Packages

library(tidyverse)
library(rstanarm)
library(bayesplot)
library(tictoc)
library(ggpubr)
library(ggrepel)
library(kableExtra)
#rstan::rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores()) # Run on multiple cores

my_ggsave <- function(plot, filename, height = 11.4, width = 11.4, units = "cm", ...){
    ggsave(filename = filename,
           plot = plot,
           height = height,
           width = width,
           units = units,
           dpi = 600,
           ...)
}

theme_set(theme_bw())

2 Data

This is how the data looks like.

d <- read_csv("data/traits.csv") %>%
  mutate(Leafhabit2 = ifelse(Leafhabit == "D", "D", "E_SD"))

d <- d %>%
  rename(log10KS = KS) %>%
  rename(log10LMA =  LMA) %>%
  rename(log10rleaf =  rleaf) %>%
  rename(log10SWCleaf =  SWCleaf) %>%
  rename(log10epsilon =  epsilon) %>%
  rename(log10Cft =  Cft) %>%
  rename(log10WOODSWC =  WOODSWC) %>%
  mutate(log10KS =  log10(log10KS)) %>%
  mutate(log10LMA =  log10(log10LMA)) %>%
  mutate(log10rleaf =  log10(log10rleaf)) %>%
  mutate(log10SWCleaf =  log10(log10SWCleaf)) %>%
  mutate(log10epsilon =  log10(log10epsilon)) %>%
  mutate(log10Cft =  log10(log10Cft)) %>%
  mutate(log10WOODSWC =  log10(log10WOODSWC))

d %>%
  head %>%
  kable() %>%
  #kable_styling("hover")
  kable_styling()
Species Sp Leafhabit Lifeform MVL log10KS P50 P88 HSM50 HSM88 HSMTLP HV rxyl log10WOODSWC po ptlp log10LMA log10rleaf Fmin log10SWCleaf log10Cft log10epsilon RWCtlp AllBranch DeadBranch Branchdiebak AllTree DeadTree Mortality Leafhabit2
Carissa spinarum CAS E S 69.83 0.1303338 -4.45 -8.59 NA NA 1.51 NA 0.72 -0.1307683 -2.38 -2.94 2.273024 -0.3467875 NA 0.0453230 -0.4317983 1.456062 88.89 NA NA NA NA NA NA E_SD
Diospyros yunnanensis DIY E T 53.40 -0.2518120 -6.41 -9.68 2.73 6.00 3.23 2.10 0.63 -0.1870866 -2.39 -3.18 2.134209 -0.1135093 -3.68 0.0211893 -0.3665315 1.158664 83.25 21 16 76.19 114 5 4.4 E_SD
Jasminum nudiflorum JAN E S 90.00 0.1367206 -5.62 -11.94 0.07 6.39 3.15 1.77 0.74 -0.2676062 -2.04 -2.47 2.025101 -0.1426675 -5.55 0.3117539 -0.2218487 1.235781 86.94 25 13 52.00 47 7 14.9 E_SD
Maytenus esquirolii MAE E S 67.83 0.0253059 -7.34 -12.34 NA NA 4.52 2.45 0.72 -0.1426675 -2.22 -2.82 2.185882 -0.2924298 -4.56 0.1430148 -0.3187588 1.304921 87.40 31 9 29.03 NA NA NA E_SD
Murraya paniculata? MUE E S 35.08 -0.2441251 -5.80 -9.35 NA NA 3.12 NA 0.83 -0.3372422 -2.06 -2.68 2.162594 NA NA 0.0086002 -0.5086383 1.442166 90.89 NA NA NA NA NA NA E_SD
Olea ferruginea OLF E T 62.50 0.0413927 -7.08 -14.15 1.45 8.52 3.57 2.49 0.87 -0.4317983 -2.98 -3.51 2.251322 -0.1135093 -5.63 0.0606978 -0.3665315 1.446692 85.12 14 0 0.00 16 0 0.0 E_SD

3 Outliers

We check outliers first.

dl <- d %>%
  pivot_longer(MVL:RWCtlp) %>%
  mutate(name2 = case_when(
  name == "log10KS" ~ "log[10]~K[S]~(kg~m^{-1}~s^{-1}~MPa^{-1})",
  name == "log10Cft" ~ "log[10]~C[leaf]~(mol~m^-2~MPa^-1)",
  name == "P50" ~ "P[50]~(MPa)",
  name == "P88" ~ "P[88]~(MPa)",
  name == "HSM50" ~ "HSM[50]~(MPa)",
  name == "HSM88" ~ "HSM[88]~(MPa)",
  name == "HSMTLP" ~ "HSM[tlp]~(MPa)",
  name == "po" ~ "pi[o]~(MPa)",
  name == "ptlp" ~ "pi[tlp]~(MPa)",
  name == "HV" ~ "H[V]~(mm^2~cm^-2)",
  name == "Fmin" ~ "Psi[min]~(MPa)",
  name == "log10epsilon" ~ "log[10]~epsilon~(MPa)",
  name == "log10LMA" ~ "log[10]~LMA~(g~m^-2)",
  name == "MVL" ~ "MVL~(cm)",
  name == "log10SWCleaf" ~ "log[10]~SWC[leaf]~(g~g^-1)",
  name == "log10WOODSWC" ~ "log[10]~SWC[wood]~(g~g^-1)",
  name == "RWCtlp" ~ "RWC[tlp]~('%')",
  name == "log10rleaf" ~ "log[10]~rho[leaf]~(g~cm^-3)",
  name == "rxyl" ~ "rho[wood]~(g~cm^-3)",
  TRUE ~ name,
                          ))


dl %>%
  filter(!is.na(DeadTree)) %>%
  ggplot(., aes(x = value, y = DeadTree/AllTree * 100, col = Leafhabit)) +
  geom_point() +
  geom_text_repel(aes(label = Sp)) +
  facet_wrap(~name, scale = "free_x") +
  ylab("Top-kill ratio (%)")

dl %>%
  filter(!is.na(DeadBranch)) %>%
  ggplot(., aes(x = value, y = DeadBranch/AllBranch * 100, col = Leafhabit)) +
  geom_point() +
  geom_text_repel(aes(label = Sp)) +
  facet_wrap(~name, scale = "free_x") +
  ylab("Branch dieback ratio (%)")

These seem to be outliers and we remove those sample from the analysis.

Trait Species
Cft WOF
epslioin PIW
HV ACFA
KS BACH
SWCleaf ARO
WOODSWC CIR
P50 PIW
P88 PIW
HSM88 PIW
HSMTLP PIW
dl <- dl %>%
  mutate(name2 = factor(name2,
         levels = c(
                    "log[10]~K[S]~(kg~m^{-1}~s^{-1}~MPa^{-1})",
                    "P[50]~(MPa)",
                    "P[88]~(MPa)",
                    "MVL~(cm)",
                    "H[V]~(mm^2~cm^-2)",
                    "Psi[min]~(MPa)",
                    "rho[wood]~(g~cm^-3)",
                    "log[10]~SWC[wood]~(g~g^-1)",
                    "pi[o]~(MPa)",
                    "pi[tlp]~(MPa)",
                    "log[10]~SWC[leaf]~(g~g^-1)",
                    "log[10]~C[leaf]~(mol~m^-2~MPa^-1)",
                    "log[10]~epsilon~(MPa)",
                    "RWC[tlp]~('%')",
                    "log[10]~LMA~(g~m^-2)",
                    "log[10]~rho[leaf]~(g~cm^-3)",
                    "HSM[50]~(MPa)",
                    "HSM[88]~(MPa)",
                    "HSM[tlp]~(MPa)"
         )))

dl0 <- dl

dl <- dl %>%
  filter(name != "log10Cft" | Sp != "WOF") %>%
  filter(name != "log10epsilon" | Sp != "PIW") %>%
  filter(name != "HV" | Sp != "ACFA") %>%
  filter(name != "log10KS" | Sp != "BACH") %>%
  filter(name != "log10SWCleaf" | Sp != "ARO") %>%
  filter(name != "log10WOODSWC" | Sp != "CIR") %>%
  filter(name != "P50" | Sp != "PIW") %>%
  filter(name != "P88" | Sp != "PIW") %>%
  filter(name != "HSM88" | Sp != "PIW") %>%
  filter(name != "HSMTLP" | Sp != "PIW")

dl_out <- dl0 %>%
  filter((name == "log10Cft" & Sp == "WOF") |
    (name == "log10epsilon" & Sp == "PIW") |
    (name == "HV" & Sp == "ACFA")  |
    (name == "log10KS" & Sp == "BACH") |
    (name == "log10SWCleaf" & Sp == "ARO")  |
    (name == "log10WOODSWC" & Sp == "CIR")  |
    (name == "P50" & Sp == "PIW")  |
    (name == "P88" & Sp == "PIW")  |
    (name == "HSM88" & Sp == "PIW")  |
    (name == "HSMTLP" & Sp == "PIW"))

4 Model

4.1 Each traits

  • We model dieback and top-kill ratios as a function of each trait.
  • We remove data that don’t have mortatliy or branch death data.
  • We also remove outliers (see above).
set.seed(123)
n_iter <- 2000
n_chain <- 4
my_prior <- normal(location = 0, scale = 10, autoscale = TRUE)
my_prior_aux <- cauchy(location = 0, scale = 2.5, autoscale = FALSE)

tic()
nd <- dl %>%
  filter(!is.na(AllBranch)) %>%
  group_by(name) %>%
  nest() %>%
#  head(2) %>%
  mutate(m_branch_LH = map(data, ~  stan_glm(cbind(DeadBranch, AllBranch - DeadBranch) ~
                                      Leafhabit*value,
                                    iter = n_iter,
                                    chains = n_chain,
                                    prior = my_prior,
                                    refresh = 0,
                                    verbose = FALSE,
                                    prior_intercept = my_prior,
                                    data = ., family = binomial))) %>%
  mutate(m_branch_LH2 = map(data, ~  stan_glm(cbind(DeadBranch, AllBranch - DeadBranch) ~
                                      Leafhabit2*value,
                                    iter = n_iter,
                                    chains = n_chain,
                                    refresh = 0,
                                    verbose = FALSE,
                                    prior = my_prior,
                                    prior_intercept = my_prior,
                                    data = ., family = binomial))) %>%
  mutate(m_branch_LF = map(data, ~  stan_glm(cbind(DeadBranch, AllBranch - DeadBranch) ~
                                      Lifeform*value,
                                    iter = n_iter,
                                    chains = n_chain,
                                    refresh = 0,
                                    verbose = FALSE,
                                    prior = my_prior,
                                    prior_intercept = my_prior,
                                    data = ., family = binomial)))

nd_tree <- dl %>%
  filter(!is.na(AllTree)) %>%
  group_by(name) %>%
  nest() %>%
  mutate(m_tree_LH = map(data, ~  stan_glm(cbind(DeadTree, AllTree - DeadTree) ~
                                      Leafhabit*value,
                                    iter = n_iter,
                                    chains = n_chain,
                                    prior = my_prior,
                                    refresh = 0,
                                    verbose = FALSE,
                                    prior_intercept = my_prior,
                                    data = ., family = binomial))) %>%
  mutate(m_tree_LH2 = map(data, ~  stan_glm(cbind(DeadTree, AllTree - DeadTree) ~
                                      Leafhabit2*value,
                                    iter = n_iter,
                                    chains = n_chain,
                                    refresh = 0,
                                    verbose = FALSE,
                                    prior = my_prior,
                                    prior_intercept = my_prior,
                                    data = ., family = binomial))) %>%
  mutate(m_tree_LF = map(data, ~  stan_glm(cbind(DeadTree, AllTree - DeadTree) ~
                                      Lifeform*value,
                                    iter = n_iter,
                                    chains = n_chain,
                                    refresh = 0,
                                    verbose = FALSE,
                                    prior = my_prior,
                                    prior_intercept = my_prior,
                                    data = ., family = binomial)))
toc()
## 231.267 sec elapsed
check_fun <- function(posterior) {
  if (quantile(posterior, 0.975) > 0 & quantile(posterior, 0.025) > 0) {
    "sig"
  } else if (quantile(posterior, 0.975) < 0 & quantile(posterior, 0.025) < 0) {
    "sig"
} else "NS"}

check_dat_fun <- function(fit, gr) {
  pred_mat <- as_tibble(fit)

  if (gr == "Leafhabit") {
  tibble(Leafhabit = c("D", "E", "SD"),
         sig_LH = c(
                 check_fun(pred_mat$value),
                 check_fun(pred_mat$value + pred_mat$`LeafhabitE:value`),
                 check_fun(pred_mat$value + pred_mat$`LeafhabitSD:value`)
                 ))
  } else if (gr == "Leafhabit2") {
  tibble(Leafhabit2 = c("D", "E_SD"),
         sig_LH2 = c(
                 check_fun(pred_mat$value),
                 check_fun(pred_mat$value + pred_mat$`Leafhabit2E_SD:value`)
                 ))
  } else if (gr == "Lifeform"){
  tibble(Lifeform = c("L", "S", "T"),
         sig_LF = c(
                 check_fun(pred_mat$value),
                 check_fun(pred_mat$value + pred_mat$`LifeformS:value`),
                 check_fun(pred_mat$value + pred_mat$`LifeformT:value`)
                 ))
  }# else {
#
#  }
}

logistic <- function(x) 1 / (1 + exp (-x))
dat_fun <- function(data, fit, gr = c("LH", "LH2", "LF")) {
  data2 <- data %>%
    dplyr::select(Leafhabit, Leafhabit2, Lifeform, value) %>%
    na.omit

  if (gr == "LH") {
    new_data <- data2 %>%
      group_by(Leafhabit) %>%
      summarize(min = min(value), max = max(value)) %>%
      group_by(Leafhabit) %>%
      nest %>%
      mutate(value = map(data, ~ seq(.[1,1] %>% unlist, .[1,2] %>% unlist, length.out = 80))) %>%
      dplyr::select(-data) %>%
      unnest(cols = c(value))
  } else if (gr == "LH2") {
    new_data <- data2 %>%
      group_by(Leafhabit2) %>%
      summarize(min = min(value), max = max(value)) %>%
      group_by(Leafhabit2) %>%
      nest %>%
      mutate(value = map(data, ~ seq(.[1,1] %>% unlist, .[1,2] %>% unlist, length.out = 80))) %>%
      dplyr::select(-data) %>%
      unnest(cols = c(value))
  } else if (gr == "LF") {
    new_data <- data2 %>%
      group_by(Lifeform) %>%
      summarize(min = min(value), max = max(value)) %>%
      group_by(Lifeform) %>%
      nest %>%
      mutate(value = map(data, ~ seq(.[1,1] %>% unlist, .[1,2] %>% unlist, length.out = 80))) %>%
      dplyr::select(-data) %>%
      unnest(cols = c(value))
  }

  pred_lin <- posterior_linpred(fit, newdata = new_data)

  tibble(
#    Leafhabit = new_data$Leafhabit,
#    Lifeform = new_data$Lifeform,
    median = pred_lin %>% apply(., 2, median),
    lower = pred_lin %>% apply(., 2, function(x)quantile(x, 0.025)),
    upper = pred_lin %>% apply(., 2, function(x)quantile(x, 0.975))) %>%
    bind_cols(new_data, .) %>%
    mutate(median_p = logistic(median)) %>%
    mutate(lower_p = logistic(lower)) %>%
    mutate(upper_p = logistic(upper))
}
nd2 <- nd %>%
  mutate(dat_branch_LH = map2(data, m_branch_LH, dat_fun, "LH")) %>%
  mutate(dat_branch_LH2 = map2(data, m_branch_LH2, dat_fun, "LH2")) %>%
  mutate(dat_branch_LF = map2(data, m_branch_LF, dat_fun, "LF"))
nd2_tree <- nd_tree %>%
  mutate(dat_tree_LH = map2(data, m_tree_LH, dat_fun, "LH")) %>%
  mutate(dat_tree_LH2 = map2(data, m_tree_LH2, dat_fun, "LH2")) %>%
  mutate(dat_tree_LF = map2(data, m_tree_LF, dat_fun, "LF"))
save.image("data/model.Rda")
load("data/model.Rda")
nd_tmp <- dl %>%
  filter(!is.na(AllBranch)) %>%
  group_by(name, name2) %>%
  nest() %>%
  dplyr::select(name, name2) %>%
  mutate(name2 = factor(name2,
         levels = c(
                    "log[10]~K[S]~(kg~m^{-1}~s^{-1}~MPa^{-1})",
                    "P[50]~(MPa)",
                    "P[88]~(MPa)",
                    "MVL~(cm)",
                    "H[V]~(mm^2~cm^-2)",
                    "Psi[min]~(MPa)",
                    "rho[wood]~(g~cm^-3)",
                    "log[10]~SWC[wood]~(g~g^-1)",
                    "pi[o]~(MPa)",
                    "pi[tlp]~(MPa)",
                    "log[10]~SWC[leaf]~(g~g^-1)",
                    "log[10]~C[leaf]~(mol~m^-2~MPa^-1)",
                    "log[10]~epsilon~(MPa)",
                    "RWC[tlp]~('%')",
                    "log[10]~LMA~(g~m^-2)",
                    "log[10]~rho[leaf]~(g~cm^-3)",
                    "HSM[50]~(MPa)",
                    "HSM[88]~(MPa)",
                    "HSM[tlp]~(MPa)"
         )
  ))

nd <- left_join(nd, nd_tmp, by = "name")
nd2 <- left_join(nd2, nd_tmp, by = "name")
nd2_tree <- left_join(nd2_tree, nd_tmp, by = "name")

4.2 Bayesian R2

Functions to get the Bayesian R2 values.

get_r2 <- function(ypred, y, median = TRUE) {
  e <- -1 * sweep(ypred, 2, y)
  var_ypred <- apply(ypred, 1, var)
  var_e <- apply(e, 1, var)
  r2 <- var_ypred / (var_ypred + var_e)
  if (median) median(r2) else r2
}

get_r2_fun <- function(fit, trait) {
  ypred <- posterior_linpred(fit) %>% logistic
  ypred <- ypred * 100
  get_r2(ypred, trait)
}

Next, we model the branch dieback and top-kill ratios as functions of groups (i.e., leaf habits or growth forms) to estimate the Bayesian R2 for the group effects alone models.

dat_branch <- d %>% filter(!is.na(DeadBranch))

b_fit_LH <- stan_glm(cbind(DeadBranch, AllBranch - DeadBranch) ~
                     Leafhabit,
                   iter = n_iter,
                   chains = n_chain,
                   prior = my_prior,
                   refresh = 0,
                   verbose = FALSE,
                   prior_intercept = my_prior,
                   data = dat_branch, family = binomial)
b_fit_LH2 <- stan_glm(cbind(DeadBranch, AllBranch - DeadBranch) ~
                     Leafhabit2,
                   iter = n_iter,
                   chains = n_chain,
                   prior = my_prior,
                   refresh = 0,
                   verbose = FALSE,
                   prior_intercept = my_prior,
                   data = dat_branch, family = binomial)
b_fit_LF <- stan_glm(cbind(DeadBranch, AllBranch - DeadBranch) ~
                     Lifeform,
                   iter = n_iter,
                   chains = n_chain,
                   prior = my_prior,
                   refresh = 0,
                   verbose = FALSE,
                   prior_intercept = my_prior,
                   data = dat_branch, family = binomial)

dat_tree <- d %>% filter(!is.na(DeadTree))

t_fit_LH <- stan_glm(cbind(DeadTree, AllTree - DeadTree) ~
                     Leafhabit,
                   iter = n_iter,
                   chains = n_chain,
                   prior = my_prior,
                   refresh = 0,
                   verbose = FALSE,
                   prior_intercept = my_prior,
                   data = dat_tree, family = binomial)
t_fit_LH2 <- stan_glm(cbind(DeadTree, AllTree - DeadTree) ~
                     Leafhabit2,
                   iter = n_iter,
                   chains = n_chain,
                   prior = my_prior,
                   refresh = 0,
                   verbose = FALSE,
                   prior_intercept = my_prior,
                   data = dat_tree, family = binomial)
t_fit_LF <- stan_glm(cbind(DeadTree, AllTree - DeadTree) ~
                     Lifeform,
                   iter = n_iter,
                   chains = n_chain,
                   prior = my_prior,
                   refresh = 0,
                   verbose = FALSE,
                   prior_intercept = my_prior,
                   data = dat_tree, family = binomial)

R2 for group effects alone models.

r2_branch_LH_v <- get_r2_fun(b_fit_LH, dat_branch$Branchdiebak)
r2_branch_LH2_v <- get_r2_fun(b_fit_LH2, dat_branch$Branchdiebak)
r2_branch_LF_v <- get_r2_fun(b_fit_LF, dat_branch$Branchdiebak)
r2_tree_LH_v <- get_r2_fun(t_fit_LH, dat_tree$Mortality)
r2_tree_LH2_v <- get_r2_fun(t_fit_LH2, dat_tree$Mortality)
r2_tree_LF_v <- get_r2_fun(t_fit_LF, dat_tree$Mortality)

r2_dat <- tibble(
       model = c(
                 "dieback-Leafhabit (E, SD, L)",
                 "dieback-Leafhabit (E + SD, L)",
                 "dieback-Growthform",
                 "top-kill-Leafhabit (E, SD, L)",
                 "top-kill-Leafhabit (E + SD, L)",
                 "top-kill-Growthform"),
       r2 = c(
              r2_branch_LH_v,
              r2_branch_LH2_v,
              r2_branch_LF_v,
              r2_tree_LH_v,
              r2_tree_LH2_v,
              r2_tree_LF_v))

r2_dat %>%
  mutate(r2 = round(r2, 2)) %>%
  kable() %>%
  kable_styling()
model r2
dieback-Leafhabit (E, SD, L) 0.28
dieback-Leafhabit (E + SD, L) 0.14
dieback-Growthform 0.05
top-kill-Leafhabit (E, SD, L) 0.34
top-kill-Leafhabit (E + SD, L) 0.20
top-kill-Growthform 0.41

For example, leaf habit alone explains 28% of variance in branch dieback (r2 = 0.28).

r2_branch <- nd2 %>%
  mutate(data2 = map(data, ~ filter(., !is.na(value))))%>%
  mutate(branch_diebak = map(data2, ~ .$Branchdiebak)) %>%
#  mutate(branch_diebak_raw = map(branch_diebak, ~ logit(. /100))) %>%
  mutate(pred_branch_LH = map(m_branch_LH, ~ logistic(posterior_linpred(.)) * 100)) %>%
  mutate(pred_branch_LH2 = map(m_branch_LH2, ~ logistic(posterior_linpred(.)) * 100)) %>%
  mutate(pred_branch_LF = map(m_branch_LF, ~ logistic(posterior_linpred(.)) * 100)) %>%
#  mutate(pred_branch_LH_raw = map(m_branch_LH, ~ posterior_linpred(.))) %>%
  mutate(r2_branch_LH = map2_dbl(pred_branch_LH, branch_diebak,  get_r2)) %>%
  mutate(r2_branch_LH2 = map2_dbl(pred_branch_LH2, branch_diebak,  get_r2)) %>%
  mutate(r2_branch_LF = map2_dbl(pred_branch_LF, branch_diebak,  get_r2)) %>%
#  mutate(r2_branch_LH_raw = map2_dbl(pred_branch_LH_raw, branch_diebak_raw,  get_r2)) %>%
  #dplyr::select(name, r2_branch_LH, r2_branch_LH_raw) %>%
  dplyr::select(name, r2_branch_LH, r2_branch_LH2, r2_branch_LF)# %>%
  #arrange(desc(r2_branch_LH))


tibble(trait = r2_branch$name, round(r2_branch[,2:4], 2)) %>%
  rename(`3 Leaf habit` = "r2_branch_LH") %>%
  rename(`2 Leaf habit` = "r2_branch_LH2") %>%
  rename(`Growth form` = "r2_branch_LF") %>%
  kable() %>%
  add_header_above(c(" ", "R2 values" = 3)) %>%
  kable_styling()
R2 values
trait 3 Leaf habit 2 Leaf habit Growth form
MVL 0.32 0.18 0.22
log10KS 0.40 0.21 0.24
P50 0.58 0.57 0.41
P88 0.68 0.64 0.34
HSM50 0.41 0.25 0.31
HSM88 0.70 0.57 0.44
HSMTLP 0.63 0.60 0.47
HV 0.25 0.16 0.15
rxyl 0.41 0.25 0.15
log10WOODSWC 0.34 0.19 0.18
po 0.50 0.46 0.41
ptlp 0.52 0.51 0.50
log10LMA 0.45 0.29 0.26
log10rleaf 0.39 0.23 0.18
Fmin 0.34 0.25 0.19
log10SWCleaf 0.45 0.34 0.22
log10Cft 0.48 0.32 0.21
log10epsilon 0.59 0.42 0.36
RWCtlp 0.47 0.36 0.34
r2_tree <- nd2_tree %>%
  mutate(data2 = map(data, ~ filter(., !is.na(value))))%>%
  mutate(mort = map(data2, ~ .$Mortality)) %>%
  mutate(pred_tree_LH = map(m_tree_LH, ~ logistic(posterior_linpred(.)) * 100)) %>%
  mutate(pred_tree_LH2 = map(m_tree_LH2, ~ logistic(posterior_linpred(.)) * 100)) %>%
  mutate(pred_tree_LF = map(m_tree_LF, ~ logistic(posterior_linpred(.)) * 100)) %>%
  mutate(r2_tree_LH = map2_dbl(pred_tree_LH, mort,  get_r2)) %>%
  mutate(r2_tree_LH2 = map2_dbl(pred_tree_LH2, mort,  get_r2)) %>%
  mutate(r2_tree_LF = map2_dbl(pred_tree_LF, mort,  get_r2)) %>%
  dplyr::select(name, r2_tree_LH, r2_tree_LH2, r2_tree_LF)# %>%
  #arrange(desc(r2_tree_LH))

tibble(trait = r2_tree$name, round(r2_tree[,2:4], 2)) %>%
  rename(`3 Leaf habit` = "r2_tree_LH") %>%
  rename(`2 Leaf habit` = "r2_tree_LH2") %>%
  rename(`Growth form` = "r2_tree_LF") %>%
  kable() %>%
  add_header_above(c(" ", "R2 values" = 3)) %>%
  kable_styling()
R2 values
trait 3 Leaf habit 2 Leaf habit Growth form
MVL 0.39 0.22 0.41
log10KS 0.42 0.18 0.38
P50 0.69 0.65 0.66
P88 0.77 0.72 0.73
HSM50 0.50 0.36 0.44
HSM88 0.65 0.55 0.57
HSMTLP 0.63 0.56 0.62
HV 0.39 0.21 0.42
rxyl 0.48 0.26 0.41
log10WOODSWC 0.35 0.19 0.64
po 0.70 0.70 0.69
ptlp 0.75 0.74 0.75
log10LMA 0.36 0.20 0.44
log10rleaf 0.75 0.39 0.74
Fmin 0.42 0.34 0.70
log10SWCleaf 0.39 0.28 0.48
log10Cft 0.74 0.63 0.61
log10epsilon 0.86 0.86 0.90
RWCtlp 0.82 0.79 0.89

5 Figures

5.1 Color code

  • Evergreen: green “#008B00”
  • Deciduous: orange “#FF8C00”
  • Semi-deciduous: blue “#1874CD”
  • Tree: light blue “#61f0f8”
  • Liana: red “#de2e19”
  • Shrub: dark blue “#1f306f”
gr <- "#008B00"
or <- "#FF8C00"
bl <- "#1874CD"
br <- "#8B7355"
lbl <- "#61F0F8"
re <- "#DE2E19"
dbl <- "#1F306F"

my_col <- c(gr, or, bl, lbl, re, dbl)

expand.grid(x = 1:3, y = 1:2) %>%
  mutate(gr = c("Evergreen", "Deciduous", "Semi-diciduous", "Tree", "Liana",
                "Shurb")) %>%
  mutate(my_col = my_col) %>%
  ggplot(aes(x, y, fill = I(my_col))) +
  geom_tile() +
  annotate("text", x = 1:3, y = 1, label = c("Evergreen", "Deciduous", "Semi-deciduous")) +
  annotate("text", x = 1:3, y = 2, label = c("Tree", "Liana", "Shrub"), size = 4) +
  coord_equal() +
  theme_void() +
  scale_y_reverse()

LF_shapes <- c("Liana"=16, "Shrub"=17, "Tree"=15)
LF_shapes2 <- c("L"=21, "S"=22, "T"=24)
LH_shapes <- c("D"=16, "E"=17, "SD"=15)

5.2 Branch dieback ratio and 3 leaf habits

sig_LH <- nd2 %>%
  mutate(sig_LH = map(m_branch_LH,  check_dat_fun, "Leafhabit")) %>%
#  mutate(sig_LF= map(m_branch_LF, check_dat_fun, "Lifeform")) %>%
  dplyr::select(sig_LH) %>%
  unnest(c( sig_LH))

sig_LH2 <- nd2 %>%
  mutate(sig_LH2 = map(m_branch_LH2,  check_dat_fun, "Leafhabit2")) %>%
#  mutate(sig_LF= map(m_branch_LF, check_dat_fun, "Lifeform")) %>%
  dplyr::select(sig_LH2) %>%
  unnest(c( sig_LH2))

sig_LF <- nd2 %>%
  mutate(sig_LF = map(m_branch_LF,  check_dat_fun, "Lifeform")) %>%
#  mutate(sig_LF= map(m_branch_LF, check_dat_fun, "Lifeform")) %>%
  dplyr::select(sig_LF) %>%
  unnest(c( sig_LF))


r2_branch2 <- r2_branch %>%
  left_join(., nd2, by = "name") %>%
  dplyr::select(name, r2_branch_LH, r2_branch_LH2, r2_branch_LF, data) %>%
  mutate(min_x = map_dbl(data, ~ min(.$value, na.rm = TRUE))) %>%
  mutate(max_x = map_dbl(data, ~ max(.$value, na.rm = TRUE))) %>%
  #dplyr::select(name, r2_branch_LH, data) %>%
  unnest(cols = data) %>%
  dplyr::select(name, name2, min_x, max_x, r2_branch_LH, r2_branch_LH2, r2_branch_LF) %>%
  unique %>%
  mutate(Leafhabit = "D") %>%
  mutate(Leafhabit2 = "D") %>%
  mutate(Lifeform = "T") %>%
  mutate(r2_branch_LF = round(r2_branch_LF, 2) %>% format(., nsmall = 2)) %>%
  mutate(r2_branch_LH2 = round(r2_branch_LH2, 2) %>% format(., nsmall = 2)) %>%
  mutate(r2_branch_LH = round(r2_branch_LH, 2) %>% format(., nsmall = 2)) %>%
  mutate(LH_lab = paste0("italic(r^2) == ", r2_branch_LH)) %>%
  mutate(LH2_lab = paste0("italic(r^2) == ", r2_branch_LH2)) %>%
  mutate(LF_lab = paste0("italic(r^2) == ", r2_branch_LF))

p_B_LH_R2 <- nd2 %>%
  dplyr::select(name2, dat_branch_LH) %>%
  unnest(c(dat_branch_LH)) %>%
  ungroup %>%
  left_join(., sig_LH, by = c("Leafhabit", "name")) %>%
  filter(sig_LH == "sig") %>%
  ggplot(., aes(x = value, col = Leafhabit, shape = Leafhabit)) +
    geom_ribbon(aes(ymin = lower_p * 100, ymax = upper_p * 100, fill =
                    Leafhabit), alpha = 0.5, col = "transparent",
                show.legend = FALSE) +
    geom_point(data = dl %>% filter(!is.na(AllBranch)), aes(x = value, y = DeadBranch/AllBranch * 100)) +
    geom_point(data = dl_out %>% filter(!is.na(AllBranch)),
               aes(y = DeadBranch/AllBranch * 100), fill = "grey", col = "black", pch = 21) +
    geom_line(aes(y = median_p * 100)) +
    geom_text(data = r2_branch2,
              aes(y = 110, x = min_x, label = LH_lab),
              hjust = -0.2,
              parse = TRUE,
              col = "black"
              ) +
    facet_wrap(~name2, scale = "free_x", label = label_parsed) +
    scale_linetype_manual(values = c(2, 1)) +
    scale_shape_manual(values = LH_shapes,
                        labels = c("Deciduous", "Evergreen", "Semi-deciduous")) +
    scale_colour_manual(values = c(or, gr, bl),
                        labels = c("Deciduous", "Evergreen", "Semi-deciduous")) +
    scale_fill_manual(values = c(or, gr, bl)) +
    scale_y_continuous(breaks = seq(0, 100, length = 5),
                       limits = c(0, 115)) +
    ylab("Branch dieback ratio (%)") +
    xlab("") +
    labs(color = "Leaf habit", shape = "Leaf habit") +
    theme(legend.position = c(0.9, 0.1))

my_ggsave(p_B_LH_R2, "figs/BD_LH_R2.png", height = 20, width = 30)
my_ggsave(p_B_LH_R2, "figs/BD_LH_R2.pdf", height = 20, width = 30)

5.3 Branch dieback ratio and 2 leaf habits

p_B_LH2_R2 <- nd2 %>%
  dplyr::select(name2, dat_branch_LH2) %>%
  unnest(c(dat_branch_LH2)) %>%
  ungroup %>%
  left_join(., sig_LH2, by = c("Leafhabit2", "name")) %>%
  filter(sig_LH2 == "sig") %>%
  ggplot(., aes(x = value, col = Leafhabit2, shape = Leafhabit2)) +
    geom_ribbon(aes(ymin = lower_p * 100, ymax = upper_p * 100, fill =
                    Leafhabit2), alpha = 0.5, col = "transparent",
                show.legend = FALSE) +
    geom_point(data = dl %>% filter(!is.na(AllBranch)), aes(x = value, y = DeadBranch/AllBranch * 100)) +
    geom_point(data = dl_out %>% filter(!is.na(AllBranch)),
               aes(y = DeadBranch/AllBranch * 100), fill = "grey", col = "black", pch = 21) +
    geom_line(aes(y = median_p * 100)) +
    geom_text(data = r2_branch2,
              aes(y = 110, x = min_x, label = LH2_lab),
              hjust = -0.2,
              parse = TRUE,
              col = "black"
              ) +
    facet_wrap(~name2, scale = "free_x", label = label_parsed) +
    scale_linetype_manual(values = c(2, 1)) +
    scale_shape_manual(values = c(16, 17),
                        labels = c("Deciduous", "Evergreen and \nSemi-deciduous")) +
    scale_colour_manual(values = c(or, gr),
                        labels = c("Deciduous", "Evergreen and \nSemi-deciduous")) +
    scale_fill_manual(values = c(or, gr)) +
    scale_y_continuous(breaks = seq(0, 100, length = 5),
                       limits = c(0, 115)) +
    ylab("Branch dieback ratio (%)") +
    xlab("") +
    labs(color = "Leaf habit", shape = "Leaf habit") +
    theme(legend.position = c(0.9, 0.1))

my_ggsave(p_B_LH2_R2, "figs/BD_LH2_R2.png", height = 20, width = 30)
my_ggsave(p_B_LH2_R2, "figs/BD_LH2_R2.pdf", height = 20, width = 30)

5.4 Branch dieback ratio and growth forms

p_B_LF_R2 <- nd2 %>%
  dplyr::select(name2, dat_branch_LF) %>%
  unnest(c(dat_branch_LF)) %>%
  ungroup %>%
  left_join(., sig_LF, by = c("Lifeform", "name")) %>%
  filter(sig_LF == "sig") %>%
  ggplot(., aes(x = value, col = Lifeform, fill = Lifeform)) +
    geom_ribbon(aes(ymin = lower_p * 100, ymax = upper_p * 100, fill =
                    Lifeform), alpha = 0.8, col = "transparent",
                show.legend = FALSE) +
    geom_point(data = dl %>% filter(!is.na(AllBranch)), aes(x = value, y = DeadBranch/AllBranch * 100, shape = Lifeform), col = "black") +
    geom_point(data = dl_out %>% filter(!is.na(AllBranch)),
               aes(y = DeadBranch/AllBranch * 100), fill = "grey", col = "black", pch = 21) +
    geom_line(aes(y = median_p * 100)) +
    geom_text(data = r2_branch2,
              aes(y = 110, x = min_x, label = LF_lab),
              hjust = -0.2,
              parse = TRUE,
              col = "black"
              ) +
    facet_wrap(~name2, scale = "free_x", label = label_parsed) +
    #scale_linetype_manual(values = c(2, 1)) +
    scale_shape_manual(values = LF_shapes2,
                      labels = c("Liana", "Shrub", "Tree")) +
    scale_colour_manual(values = c(re, dbl, lbl
                                #)) +
                                ), guide = "none") +
    scale_fill_manual(values = c(re, dbl, lbl),
                      labels = c("Liana", "Shrub", "Tree")
                                ) +
    ylab("Branch dieback ratio (%)") +
    xlab("") +
    labs(fill = "Growth form", shape = "Growth form", col = "Growth form") +
    scale_y_continuous(breaks = seq(0, 100, length = 5),
                       limits = c(0, 115)) +
    theme(legend.position = c(0.9, 0.1))

my_ggsave(p_B_LF_R2, "figs/BD_LF_R2.png", height = 20, width = 30)
my_ggsave(p_B_LF_R2, "figs/BD_LF_R2.pdf", height = 20, width = 30)

5.5 Top-kill ratio and 3 leaf habits

sig_LH <- nd2_tree %>%
  mutate(sig_LH = map(m_tree_LH,  check_dat_fun, "Leafhabit")) %>%
#  mutate(sig_LF= map(m_tree_LF, check_dat_fun, "Lifeform")) %>%
  dplyr::select(sig_LH) %>%
  unnest(c( sig_LH))

sig_LH2 <- nd2_tree %>%
  mutate(sig_LH2 = map(m_tree_LH2,  check_dat_fun, "Leafhabit2")) %>%
#  mutate(sig_LF= map(m_tree_LF, check_dat_fun, "Lifeform")) %>%
  dplyr::select(sig_LH2) %>%
  unnest(c( sig_LH2))

sig_LF <- nd2_tree %>%
  mutate(sig_LF = map(m_tree_LF,  check_dat_fun, "Lifeform")) %>%
#  mutate(sig_LF= map(m_tree_LF, check_dat_fun, "Lifeform")) %>%
  dplyr::select(sig_LF) %>%
  unnest(c( sig_LF))

r2_tree2 <- r2_tree %>%
  left_join(., nd2, by = "name") %>%
  dplyr::select(name, r2_tree_LH, r2_tree_LH2, r2_tree_LF, data) %>%
  mutate(min_x = map_dbl(data, ~ min(.$value, na.rm = TRUE))) %>%
  mutate(max_x = map_dbl(data, ~ max(.$value, na.rm = TRUE))) %>%
  #dplyr::select(name, r2_tree_LH, data) %>%
  unnest(cols = data) %>%
  dplyr::select(name, name2, min_x, max_x, r2_tree_LH, r2_tree_LH2, r2_tree_LF) %>%
  unique %>%
  mutate(Leafhabit = "D") %>%
  mutate(Leafhabit2 = "D") %>%
  mutate(Lifeform = "T") %>%
  mutate(r2_tree_LF = round(r2_tree_LF, 2) %>% format(., nsmall = 2)) %>%
  mutate(r2_tree_LH2 = round(r2_tree_LH2, 2) %>% format(., nsmall = 2)) %>%
  mutate(r2_tree_LH = round(r2_tree_LH, 2) %>% format(., nsmall = 2)) %>%
  mutate(LH_lab = paste0("italic(r^2) == ", r2_tree_LH)) %>%
  mutate(LH2_lab = paste0("italic(r^2) == ", r2_tree_LH2)) %>%
  mutate(LF_lab = paste0("italic(r^2) == ", r2_tree_LF))

p_T_LH_R2 <- nd2_tree %>%
  dplyr::select(name2, dat_tree_LH) %>%
  unnest(c(dat_tree_LH)) %>%
  ungroup %>%
  left_join(., sig_LH, by = c("Leafhabit", "name")) %>%
  filter(sig_LH == "sig") %>%
  ggplot(., aes(x = value, col = Leafhabit, shape = Leafhabit)) +
    geom_ribbon(aes(ymin = lower_p * 100, ymax = upper_p * 100, fill =
                    Leafhabit), alpha = 0.5, col = "transparent",
                show.legend = FALSE) +
    geom_point(data = dl %>% filter(!is.na(AllTree)), aes(x = value, y = DeadTree/AllTree * 100)) +
    geom_point(data = dl_out %>% filter(!is.na(AllTree)),
               aes(y = DeadTree/AllTree * 100), fill = "grey", col = "black", pch = 21) +
    geom_line(aes(y = median_p * 100)) +
    geom_text(data = r2_tree2,
              aes(y = 110, x = min_x, label = LH_lab),
              hjust = -0.2,
              parse = TRUE,
              col = "black"
              ) +
    facet_wrap(~name2, scale = "free_x", label = label_parsed) +
    scale_linetype_manual(values = c(2, 1)) +
    scale_shape_manual(values = LH_shapes,
                        labels = c("Deciduous", "Evergreen", "Semi-deciduous")) +
    scale_colour_manual(values = c(or, gr, bl),
                        labels = c("Deciduous", "Evergreen", "Semi-deciduous")) +
    scale_fill_manual(values = c(or, gr, bl)) +
    scale_y_continuous(breaks = seq(0, 100, length = 5),
                       limits = c(0, 115)) +
    ylab("Top-kill ratio (%)") +
    xlab("") +
    labs(color = "Leaf habit", shape = "Leaf habit") +
    theme(legend.position = c(0.9, 0.1))

my_ggsave(p_T_LH_R2, "figs/TD_LH_R2.png", height = 20, width = 30)
my_ggsave(p_T_LH_R2, "figs/TD_LH_R2.pdf", height = 20, width = 30)

5.6 Top-kill ratio and 2 leaf habits

p_T_LH2_R2 <- nd2_tree %>%
  dplyr::select(name2, dat_tree_LH2) %>%
  unnest(c(dat_tree_LH2)) %>%
  ungroup %>%
  left_join(., sig_LH2, by = c("Leafhabit2", "name")) %>%
  filter(sig_LH2 == "sig") %>%
  ggplot(., aes(x = value, col = Leafhabit2, shape = Leafhabit2)) +
    geom_ribbon(aes(ymin = lower_p * 100, ymax = upper_p * 100, fill =
                    Leafhabit2), alpha = 0.5, col = "transparent", show.legend
  = FALSE) +
    geom_point(data = dl %>% filter(!is.na(AllTree)), aes(x = value, y = DeadTree/AllTree * 100)) +
    geom_point(data = dl_out %>% filter(!is.na(AllTree)),
               aes(y = DeadTree/AllTree * 100), fill = "grey", col = "black", pch = 21) +
    geom_line(aes(y = median_p * 100)) +
    geom_text(data = r2_tree2,
              aes(y = 110, x = min_x, label = LH_lab),
              hjust = -0.2,
              parse = TRUE,
              col = "black"
              ) +
    facet_wrap(~name2, scale = "free_x", label = label_parsed) +
    scale_linetype_manual(values = c(2, 1)) +
    scale_shape_manual(values = c(16, 17),
                        labels = c("Deciduous", "Evergreen and \nSemi-deciduous")) +
    scale_colour_manual(values = c(or, gr),
                        labels = c("Deciduous", "Evergreen and \nSemi-deciduous")) +
    scale_fill_manual(values = c(or, gr
                                )) +
    scale_y_continuous(breaks = seq(0, 100, length = 5),
                       limits = c(0, 115)) +
    ylab("Top-kill ratio (%)") +
    xlab("") +
    labs(color = "Leaf habit", shape = "Leaf habit") +
    theme(legend.position = c(0.9, 0.1))

my_ggsave(p_T_LH2_R2, "figs/TD_LH2_R2.png", height = 20, width = 25)
my_ggsave(p_T_LH2_R2, "figs/TD_LH2_R2.pdf", height = 20, width = 25)

5.7 Top-kill ratio and growth forms

p_T_LF_R2 <- nd2_tree %>%
  dplyr::select(name2, dat_tree_LF) %>%
  unnest(c(dat_tree_LF)) %>%
  ungroup %>%
  left_join(., sig_LF, by = c("Lifeform", "name")) %>%
  filter(sig_LF == "sig") %>%
  ggplot(., aes(x = value, col = Lifeform, fill = Lifeform)) +
    geom_ribbon(aes(ymin = lower_p * 100, ymax = upper_p * 100, fill =
                    Lifeform), alpha = 0.8, col = "transparent",
                show.legend = FALSE) +
    geom_point(data = dl %>% filter(!is.na(AllTree)),
               aes(x = value, y = DeadTree/AllTree * 100,
                   shape = Lifeform),
               col = "black"
               ) +
    geom_point(data = dl_out %>% filter(!is.na(AllTree)),
               aes(y = DeadTree/AllTree * 100), fill = "grey", col = "black", pch = 21) +
    geom_line(aes(y = median_p * 100)) +
    geom_text(data = r2_tree2,
              aes(y = 110, x = min_x, label = LH_lab),
              hjust = -0.2,
              parse = TRUE,
              col = "black"
              ) +
    facet_wrap(~name2, scale = "free_x", label = label_parsed) +
    scale_shape_manual(values = LF_shapes2,
                      labels = c("Liana", "Shrub", "Tree")) +
    scale_colour_manual(values = c(re, dbl, lbl
                                ), guide = "none") +
    scale_fill_manual(values = c(re, dbl, lbl),
                      labels = c("Liana", "Shrub", "Tree")
                                ) +
    labs(fill = "Growth form") +
    scale_y_continuous(breaks = seq(0, 100, length = 5),
                       limits = c(0, 115)) +
    labs(fill = "Growth form", shape = "Growth form", col = "Growth form") +
    ylab("Top-kill ratio (%)") +
    xlab("") +
    theme(legend.position = c(0.9, 0.1))

my_ggsave(p_T_LF_R2, "figs/TD_LF_R2.png", height = 20, width = 30)
my_ggsave(p_T_LF_R2, "figs/TD_LF_R2.pdf", height = 20, width = 30)

6 Computing Environment

devtools::session_info()
## ─ Session info ───────────────────────────────────────────────────────────────
##  setting  value                       
##  version  R version 4.0.5 (2021-03-31)
##  os       Ubuntu 20.04.2 LTS          
##  system   x86_64, linux-gnu           
##  ui       X11                         
##  language (EN)                        
##  collate  en_US.UTF-8                 
##  ctype    en_US.UTF-8                 
##  tz       Etc/UTC                     
##  date     2021-07-15                  
## 
## ─ Packages ───────────────────────────────────────────────────────────────────
##  package      * version  date       lib source        
##  abind          1.4-5    2016-07-21 [1] RSPM (R 4.0.3)
##  assertthat     0.2.1    2019-03-21 [1] RSPM (R 4.0.3)
##  backports      1.2.1    2020-12-09 [1] RSPM (R 4.0.3)
##  base64enc      0.1-3    2015-07-28 [1] RSPM (R 4.0.3)
##  bayesplot    * 1.8.0    2021-01-10 [1] RSPM (R 4.0.3)
##  boot           1.3-27   2021-02-12 [2] CRAN (R 4.0.5)
##  broom          0.7.6    2021-04-05 [1] RSPM (R 4.0.4)
##  bslib          0.2.5.1  2021-05-18 [1] RSPM (R 4.0.5)
##  cachem         1.0.5    2021-05-15 [1] RSPM (R 4.0.4)
##  callr          3.6.0    2021-03-28 [1] RSPM (R 4.0.4)
##  car            3.0-10   2020-09-29 [1] RSPM (R 4.0.3)
##  carData        3.0-4    2020-05-22 [1] RSPM (R 4.0.3)
##  cellranger     1.1.0    2016-07-27 [1] RSPM (R 4.0.3)
##  cli            3.0.0    2021-06-30 [1] RSPM (R 4.0.5)
##  codetools      0.2-18   2020-11-04 [2] CRAN (R 4.0.5)
##  colorspace     2.0-0    2020-11-11 [1] RSPM (R 4.0.3)
##  colourpicker   1.1.0    2020-09-14 [1] RSPM (R 4.0.2)
##  crayon         1.4.1    2021-02-08 [1] RSPM (R 4.0.3)
##  crosstalk      1.1.1    2021-01-12 [1] RSPM (R 4.0.3)
##  curl           4.3      2019-12-02 [1] RSPM (R 4.0.3)
##  data.table     1.14.0   2021-02-21 [1] RSPM (R 4.0.3)
##  DBI            1.1.1    2021-01-15 [1] RSPM (R 4.0.3)
##  dbplyr         2.1.1    2021-04-06 [1] RSPM (R 4.0.4)
##  desc           1.3.0    2021-03-05 [1] RSPM (R 4.0.3)
##  devtools       2.4.0    2021-04-07 [1] RSPM (R 4.0.4)
##  digest         0.6.27   2020-10-24 [1] RSPM (R 4.0.3)
##  dplyr        * 1.0.7    2021-06-18 [1] RSPM (R 4.0.5)
##  DT             0.18     2021-04-14 [1] RSPM (R 4.0.4)
##  dygraphs       1.1.1.6  2018-07-11 [1] RSPM (R 4.0.3)
##  ellipsis       0.3.2    2021-04-29 [1] RSPM (R 4.0.4)
##  evaluate       0.14     2019-05-28 [1] RSPM (R 4.0.3)
##  fansi          0.5.0    2021-05-25 [1] RSPM (R 4.0.5)
##  farver         2.1.0    2021-02-28 [1] RSPM (R 4.0.3)
##  fastmap        1.1.0    2021-01-25 [1] RSPM (R 4.0.3)
##  forcats      * 0.5.1    2021-01-27 [1] RSPM (R 4.0.3)
##  foreign        0.8-81   2020-12-22 [2] CRAN (R 4.0.5)
##  fs             1.5.0    2020-07-31 [1] RSPM (R 4.0.3)
##  generics       0.1.0    2020-10-31 [1] RSPM (R 4.0.3)
##  ggplot2      * 3.3.3    2020-12-30 [1] RSPM (R 4.0.3)
##  ggpubr       * 0.4.0    2020-06-27 [1] RSPM (R 4.0.3)
##  ggrepel      * 0.9.1    2021-01-15 [1] RSPM (R 4.0.3)
##  ggridges       0.5.3    2021-01-08 [1] RSPM (R 4.0.3)
##  ggsignif       0.6.1    2021-02-23 [1] RSPM (R 4.0.3)
##  glue           1.4.2    2020-08-27 [1] RSPM (R 4.0.3)
##  gridExtra      2.3      2017-09-09 [1] RSPM (R 4.0.3)
##  gtable         0.3.0    2019-03-25 [1] RSPM (R 4.0.3)
##  gtools         3.8.2    2020-03-31 [1] RSPM (R 4.0.3)
##  haven          2.4.1    2021-04-23 [1] RSPM (R 4.0.5)
##  highr          0.9      2021-04-16 [1] RSPM (R 4.0.4)
##  hms            1.1.0    2021-05-17 [1] RSPM (R 4.0.4)
##  htmltools      0.5.1.1  2021-01-22 [1] RSPM (R 4.0.3)
##  htmlwidgets    1.5.3    2020-12-10 [1] RSPM (R 4.0.3)
##  httpuv         1.6.1    2021-05-07 [1] RSPM (R 4.0.5)
##  httr           1.4.2    2020-07-20 [1] RSPM (R 4.0.3)
##  igraph         1.2.6    2020-10-06 [1] RSPM (R 4.0.3)
##  inline         0.3.17   2020-12-01 [1] RSPM (R 4.0.3)
##  jquerylib      0.1.4    2021-04-26 [1] RSPM (R 4.0.4)
##  jsonlite       1.7.2    2020-12-09 [1] RSPM (R 4.0.3)
##  kableExtra   * 1.3.4    2021-02-20 [1] RSPM (R 4.0.3)
##  knitr          1.33     2021-04-24 [1] RSPM (R 4.0.4)
##  labeling       0.4.2    2020-10-20 [1] RSPM (R 4.0.3)
##  later          1.2.0    2021-04-23 [1] RSPM (R 4.0.4)
##  lattice        0.20-41  2020-04-02 [2] CRAN (R 4.0.5)
##  lifecycle      1.0.0    2021-02-15 [1] RSPM (R 4.0.3)
##  lme4           1.1-26   2020-12-01 [1] RSPM (R 4.0.3)
##  loo            2.4.1    2020-12-09 [1] RSPM (R 4.0.3)
##  lubridate      1.7.10   2021-02-26 [1] RSPM (R 4.0.3)
##  magrittr       2.0.1    2020-11-17 [1] RSPM (R 4.0.3)
##  markdown       1.1      2019-08-07 [1] RSPM (R 4.0.3)
##  MASS           7.3-53.1 2021-02-12 [2] CRAN (R 4.0.5)
##  Matrix         1.3-2    2021-01-06 [2] CRAN (R 4.0.5)
##  matrixStats    0.58.0   2021-01-29 [1] RSPM (R 4.0.3)
##  memoise        2.0.0    2021-01-26 [1] RSPM (R 4.0.3)
##  mime           0.11     2021-06-23 [1] RSPM (R 4.0.5)
##  miniUI         0.1.1.1  2018-05-18 [1] RSPM (R 4.0.3)
##  minqa          1.2.4    2014-10-09 [1] RSPM (R 4.0.3)
##  modelr         0.1.8    2020-05-19 [1] RSPM (R 4.0.3)
##  munsell        0.5.0    2018-06-12 [1] RSPM (R 4.0.3)
##  nlme           3.1-152  2021-02-04 [2] CRAN (R 4.0.5)
##  nloptr         1.2.2.2  2020-07-02 [1] RSPM (R 4.0.3)
##  openxlsx       4.2.3    2020-10-27 [1] RSPM (R 4.0.3)
##  pillar         1.6.1    2021-05-16 [1] RSPM (R 4.0.4)
##  pkgbuild       1.2.0    2020-12-15 [1] RSPM (R 4.0.3)
##  pkgconfig      2.0.3    2019-09-22 [1] RSPM (R 4.0.3)
##  pkgload        1.2.1    2021-04-06 [1] RSPM (R 4.0.4)
##  plyr           1.8.6    2020-03-03 [1] RSPM (R 4.0.3)
##  prettyunits    1.1.1    2020-01-24 [1] RSPM (R 4.0.3)
##  processx       3.5.1    2021-04-04 [1] RSPM (R 4.0.4)
##  promises       1.2.0.1  2021-02-11 [1] RSPM (R 4.0.3)
##  ps             1.6.0    2021-02-28 [1] RSPM (R 4.0.3)
##  purrr        * 0.3.4    2020-04-17 [1] RSPM (R 4.0.3)
##  R6             2.5.0    2020-10-28 [1] RSPM (R 4.0.3)
##  Rcpp         * 1.0.7    2021-07-07 [1] RSPM (R 4.0.5)
##  RcppParallel   5.1.2    2021-04-15 [1] RSPM (R 4.0.4)
##  readr        * 1.4.0    2020-10-05 [1] RSPM (R 4.0.4)
##  readxl         1.3.1    2019-03-13 [1] RSPM (R 4.0.3)
##  remotes        2.3.0    2021-04-01 [1] RSPM (R 4.0.4)
##  reprex         2.0.0    2021-04-02 [1] RSPM (R 4.0.4)
##  reshape2       1.4.4    2020-04-09 [1] RSPM (R 4.0.3)
##  rio            0.5.26   2021-03-01 [1] RSPM (R 4.0.3)
##  rlang          0.4.11   2021-04-30 [1] RSPM (R 4.0.4)
##  rmarkdown      2.7      2021-02-19 [1] RSPM (R 4.0.3)
##  rprojroot      2.0.2    2020-11-15 [1] RSPM (R 4.0.3)
##  rsconnect      0.8.17   2021-04-09 [1] RSPM (R 4.0.4)
##  rstan          2.21.2   2020-07-27 [1] RSPM (R 4.0.4)
##  rstanarm     * 2.21.1   2020-07-20 [1] RSPM (R 4.0.4)
##  rstantools     2.1.1    2020-07-06 [1] RSPM (R 4.0.3)
##  rstatix        0.7.0    2021-02-13 [1] RSPM (R 4.0.3)
##  rstudioapi     0.13     2020-11-12 [1] RSPM (R 4.0.3)
##  rvest          1.0.0    2021-03-09 [1] RSPM (R 4.0.3)
##  sass           0.4.0    2021-05-12 [1] RSPM (R 4.0.4)
##  scales         1.1.1    2020-05-11 [1] RSPM (R 4.0.3)
##  sessioninfo    1.1.1    2018-11-05 [1] RSPM (R 4.0.3)
##  shiny          1.6.0    2021-01-25 [1] RSPM (R 4.0.3)
##  shinyjs        2.0.0    2020-09-09 [1] RSPM (R 4.0.3)
##  shinystan      2.5.0    2018-05-01 [1] RSPM (R 4.0.0)
##  shinythemes    1.2.0    2021-01-25 [1] RSPM (R 4.0.3)
##  StanHeaders    2.21.0-7 2020-12-17 [1] RSPM (R 4.0.4)
##  statmod        1.4.35   2020-10-19 [1] RSPM (R 4.0.3)
##  stringi        1.6.2    2021-05-17 [1] RSPM (R 4.0.4)
##  stringr      * 1.4.0    2019-02-10 [1] RSPM (R 4.0.3)
##  survival       3.2-10   2021-03-16 [2] CRAN (R 4.0.5)
##  svglite        2.0.0    2021-02-20 [1] RSPM (R 4.0.4)
##  systemfonts    1.0.1    2021-02-09 [1] RSPM (R 4.0.4)
##  testthat       3.0.2    2021-02-14 [1] RSPM (R 4.0.3)
##  threejs        0.3.3    2020-01-21 [1] RSPM (R 4.0.0)
##  tibble       * 3.1.2    2021-05-16 [1] RSPM (R 4.0.4)
##  tictoc       * 1.0.1    2021-04-19 [1] RSPM (R 4.0.4)
##  tidyr        * 1.1.3    2021-03-03 [1] RSPM (R 4.0.4)
##  tidyselect     1.1.1    2021-04-30 [1] RSPM (R 4.0.4)
##  tidyverse    * 1.3.1    2021-04-15 [1] RSPM (R 4.0.4)
##  usethis        2.0.1    2021-02-10 [1] RSPM (R 4.0.3)
##  utf8           1.2.1    2021-03-12 [1] RSPM (R 4.0.3)
##  V8             3.4.1    2021-04-23 [1] RSPM (R 4.0.4)
##  vctrs          0.3.8    2021-04-29 [1] RSPM (R 4.0.4)
##  viridisLite    0.4.0    2021-04-13 [1] RSPM (R 4.0.4)
##  webshot        0.5.2    2019-11-22 [1] RSPM (R 4.0.3)
##  withr          2.4.2    2021-04-18 [1] RSPM (R 4.0.4)
##  xfun           0.24     2021-06-15 [1] RSPM (R 4.0.5)
##  xml2           1.3.2    2020-04-23 [1] RSPM (R 4.0.3)
##  xtable         1.8-4    2019-04-21 [1] RSPM (R 4.0.3)
##  xts            0.12.1   2020-09-09 [1] RSPM (R 4.0.3)
##  yaml           2.2.1    2020-02-01 [1] RSPM (R 4.0.3)
##  zip            2.1.1    2020-08-27 [1] RSPM (R 4.0.3)
##  zoo            1.8-9    2021-03-09 [1] RSPM (R 4.0.3)
## 
## [1] /usr/local/lib/R/site-library
## [2] /usr/local/lib/R/library
## [3] /home/mattocci/R/x86_64-pc-linux-gnu-library/4.0